/**                   GNU GENERAL PUBLIC LICENSE
 *                       Version 3, 29 June 2007
 *
 * Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>
 * Everyone is permitted to copy and distribute verbatim copies
 * of this license document, but changing it is not allowed.
 * See License.txt file that comes with this distribution
 *
 *
 * Author: Vitaly German <vitaly.german@gmail.com>, (C) 2007
 * 
 **/

#pragma once
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/lu.hpp>

/// Solves system of equlities using lu decomposition method
/// The input system may contain redundant inequlities
/**
 * colluns in the input matrix corresponds to the output vector
 * MX(equalities, variables)
 * C(constants)
 * Invariant: MX.size1() == C.size()
 * Algorithm:
 * 
 */
/*
template <class MX, class C>
class  EqSimpl {
    MX & mx;
    C & cnt;
public:
    
    EqSolve(MX& m, C & c):mx(m),cnt(c){}
    
    void Accept(MX& mx, C & c);
    void Solve(MX& m, C & c, size_t line_id, size_t var_id);
    void Solve() ;
    
    typedef typename C::value_type value_type;
private:
    
    typedef std::pair<size_t, value_type> Pos2Coef;
    typedef std::vector<Pos2Coef> VecPos2Coef;
    
    VecPos2Coef& build_signs(VecPos2Coef&,const size_t& line_id, const size_t& var_id) const;
    
};
*/
/**
 * Solves linear equtions for arbitrary matrix
 * Ax=b
 * a - input matrix A
 * b - input vector of constants b
 * output:
 *      true and vector of solutions x. Solutions are placed in the vector b.
 *      false if the system has no solutions or redundant
 *
 **/
template < typename M, typename V>
bool solve_equ(M & a, V & b ){
	using namespace boost::numeric::ublas;
    typedef permutation_matrix<std::size_t> pmatrix;
    pmatrix P(a.size1());
    int ret = lu_factorize(a,P);
    
    if (ret != 0){
    	b.clear();
    	return false;
    }
    
    lu_substitute(a,P,b);
    return true;
}
